


##### Packages ######
library(ggborderline)
library(shadowtext)
library(tidyverse)
library(miceadds)
library(interplot)
library(estimatr)
library(ggpubr)
library(haven)
library(countrycode)
library(plyr)
library(tidyverse)
library(readxl)
library(sjlabelled)
library(modelsummary)
library(stdmod)
library(scales)
library(gt)
library(PanelMatch)

# Theme
theme_r <- function() {
  back="white"
  font_size=9.5
  theme_minimal() +
    theme(text=element_text(family="Helvetica",size=font_size,color="black"),
          axis.text=element_text(color="black",size=font_size),
          axis.title = element_text(size=font_size),
          axis.title.x = element_text(margin = margin(3,0,0,0)),
          plot.margin = margin(10,10,10,10),
          plot.title = element_text(family="Helvetica",
                                    hjust=0,
                                    face = "bold",
                                    size=font_size),
          plot.subtitle = element_text(family="Helvetica",
                                       lineheight = 1,
                                       size=font_size,
                                       face = "italic",
                                       hjust=0,
                                       margin = margin(0,0,10,0)),
          panel.background = element_rect(fill=back,
                                          color=back),
          plot.background = element_rect(fill=back,
                                         color=back),
          panel.grid.major = element_line(color="grey95"),
          panel.grid = element_blank(),
          axis.ticks = element_line(),
          axis.line = element_line(),
          legend.background = element_rect(fill=back),
          plot.title.position = "plot",
          strip.placement = "outside",
          strip.text = element_text(size=font_size),
          plot.caption = element_text(hjust=0,size=font_size,
                                      color="grey40",
                                      lineheight = 1),
          legend.text = element_text(size=font_size),
          plot.caption.position = "plot",
          update_geom_defaults("borderline", list(size = 3)),
          update_geom_defaults("path", list(size = 3)),
          update_geom_defaults("point", list(size = 3)),
          update_geom_defaults("shadowtext", list(lineheight = 1,
                                                  bg.color=back,
                                                  color="black",
                                                  family="Helvetica")),
          update_geom_defaults("text", list(lineheight = 1,
                                            size=font_size/(14/5),
                                            family="Helvetica")),
          update_geom_defaults("line", list(linewidth = 1)))
} 
library(ggsci)
pal <- pal_npg()(10)



######################################################
############   Revise and Resubmit R3     ############ 
######################################################

prr<-data.frame(read_excel("country-year data.xlsx"))

# Construct additional variables
prr <- prr %>%
  group_by(code) %>%
  arrange(year) %>%
  dplyr::mutate(PRR_government_support_lag1 = dplyr::lag(PRR_government_support,1),
                prr_in_government_total_years = sum(PRR_government_support,na.rm=T),
                PRR_parl_dummy_lag1 = dplyr::lag(PRR_parl_dummy,1),
                prr_in_parl_total_years = sum(PRR_parl_dummy,na.rm=T))

## Vote change from previous to most recent election
library(purrr)

# Lag: previous unique value 
lag_diff <- function(x) {
  r <- rle(x)
  rep(c(NA, head(r$values, -1)), r$lengths)
}

lag_diff_0 <- function(x) {
  r <- rle(x)
  lag_values <- c(NA, head(r$values, -1))
  result <- rep(lag_values, r$lengths)
  # replace NA where no prior > 0 exists
  if (all(x <= 0)) {
    result[] <- 0
  } else {
    first_positive_idx <- which(x > 0)[1]
    result[seq_len(first_positive_idx)] <- 0
  }
  result
}

prr <- prr %>%
  arrange(country, year) %>%
  group_by(country) %>%
  dplyr::mutate(
    PRR_vote_share_prev_elect = ifelse(PRR_vote_share==0,
                                       lag_diff_0(PRR_vote_share),
                                       lag_diff(PRR_vote_share)),
    PRR_vote_share_growth = PRR_vote_share-PRR_vote_share_prev_elect,
  ) %>%
  ungroup()

# Claasen data on satisfaction with democracy
cl <- read_csv("satis_est_v2.csv")
cl$code <- cl$ISO_code
cl$year <- cl$Year

# Combine with PRR data
cl <- merge(prr,cl,by=c("code","year"),all.x = T)
cl$country[cl$country=="United Kingdom"] <- "UK"
cl$country[cl$country=="Czech Republic"] <- "Czechia"
cl <- cl %>% filter(!is.na(Satis))


################# Treatment History Figure #################
f1<-ggplot(cl,aes(year,Satis)) +
  geom_vline(aes(xintercept=year,
                 color=factor(PRR_government_support)),
             size=0.5) +  
  geom_line() +
  theme_r() +
  facet_wrap(~reorder(country,-PRR_government_support)) +
  theme(legend.position = "none",
        axis.text.x = element_text(angle=30,hjust=1)) +
  scale_color_manual(values=c("transparent","salmon")) +
  scale_x_continuous(breaks=c(1980,2000,2020)) +
  labs(x=NULL,
       y="Satisfaction with democracy")
ggsave(f1,
       file="PRR, treatment history.jpeg",
       width=6,height=6,
       dpi=1000)



################# Regression analysis #################

# TABLE: REGRESSION TABLE
m1 <- lm_robust(Satis~PRR_government_support_lag1+
                  code,
                data=cl,
                clusters = code,
                se_type="stata")
m2 <- lm_robust(Satis~PRR_government_support_lag1 +
                  vdem_corr + 
                  unemployment +
                  wdi_gdpgr + 
                  wpp_netmig+
                  code,
                data=cl,
                clusters = code,
                se_type="stata")
m3 <- lm_robust(Satis~PRR_parl_dummy_lag1+
                  code,
                data=cl,
                clusters = code,
                se_type="stata")
m4 <- lm_robust(Satis~PRR_parl_dummy_lag1 +
                  vdem_corr + 
                  unemployment +
                  wdi_gdpgr + 
                  wpp_netmig+
                  code,
                data=cl,
                clusters = code,
                se_type="stata")

modelsummary(list(m1,m2,m3,m4),
             output = "gt",
             stars = T, 
             coef_map = c("PRR_government_support_lag1"="PRR government inclusion",
                          "PRR_parl_dummy_lag1" = "PRR parliament inclusion",
                          "vdem_corr" = "Corruption (0-1)",
                          "unemployment" = "Unemployment (%)",
                          "wdi_gdpgr" = "GDP growth (%)",
                          "wpp_netmig" = "Net immigration (per 1000)",
                          "(Intercept)" = "Intercept"),
             gof_map = c("nobs","r.squared")) %>%
  tab_spanner(label="Dependent variable: Satisfaction with democracy",
              columns=c(2:5)) %>%  
  gtsave(filename="prr reg table.docx")




########## Panel matching analysis ########## 

# GOVERNMENT
cl$id <- as.integer(as.numeric(as.factor(cl$country)))
cl$time <- as.integer(cl$year)
method <- "ps.weight"
## Government support start
cl1 <- PanelData(panel.data = cl,
                 unit.id = "id",
                 time.id = "time",
                 treatment = "PRR_government_support",
                 outcome = "Satis")
p_govsup_att <- PanelMatch(cl1,
                           lag = 5,
                           refinement.method = method, 
                           match.missing = TRUE,
                           covs.formula = ~ vdem_corr + wdi_gdpgr + unemployment + wpp_netmig + I(lag(Satis, 1:4)),
                           size.match = 5, qoi = "att", 
                           lead = 0:10, forbid.treatment.reversal = FALSE,
                           use.diagonal.variance.matrix = TRUE)
dat_att<-PanelEstimate(cl1, sets = p_govsup_att,
                       number.iterations = 5000,
                       confidence.level = .95)
dat1<-data.frame(est=dat_att$estimate,
                 se=dat_att$standard.error)
dat1$time <- rownames(dat1)
dat1$conf.hi <- dat1$est+dat1$se*1.96
dat1$conf.lo <- dat1$est-dat1$se*1.96
dat1$order <- seq(0,10,1)
f1_gov<-ggplot(dat1,aes(reorder(time,order),est,group=1)) +
  geom_hline(yintercept = 0,linetype="dashed") +
  # geom_line() +
  geom_point() +
  geom_errorbar(aes(ymin=conf.lo,
                    ymax=conf.hi),width=0) +
  theme_r() +
  scale_y_continuous(limits=c(-0.5,1)) +  
  scale_x_discrete(breaks=c("t+0","t+2","t+4","t+6","t+8","t+10")) +
  labs(y="Average treatment effect",
       x=NULL,
       title="Government entry")
### Tables for matched sets
cnames <- cl %>% 
  group_by(country) %>%
  dplyr::select(id,country) %>%
  slice(1)
### Entry
att_match<-list(dat_att$matched.sets)
att_match<-data.frame(id_year = objects(att_match[[1]]))
att_match<-att_match %>% mutate(id=gsub("\\..*","",id_year),
                                year=gsub("^.*\\.","",id_year))
att_match <- merge(att_match,cnames,by="id",all.x=T)
att_match <- att_match %>%
  mutate(cy = paste0(country," (",year,")")) %>%
  dplyr::select(cy)
paste(att_match$cy,collapse = " ")

## Government support end
cl1 <- PanelData(panel.data = cl,
                 unit.id = "id",
                 time.id = "time",
                 treatment = "PRR_government_support",
                 outcome = "Satis")
p_govsup_art <- PanelMatch(cl1,
                           lag = 4,
                           refinement.method = method, 
                           match.missing = TRUE,
                           covs.formula = ~ vdem_corr + wdi_gdpgr + unemployment + wpp_netmig + I(lag(Satis, 1:2)),
                           size.match = 5, qoi = "art",
                           lead = 0:5, forbid.treatment.reversal = FALSE,
                           use.diagonal.variance.matrix = TRUE)
dat_art<-PanelEstimate(cl1, sets = p_govsup_art,
                       number.iterations = 5000,
                       confidence.level = .95)
dat1<-data.frame(est=dat_art$estimate,
                 se=dat_art$standard.error)
dat1$time <- rownames(dat1)
dat1$conf.hi <- dat1$est+dat1$se*1.96
dat1$conf.lo <- dat1$est-dat1$se*1.96
dat1$order <- seq(0,5,1)
f2_gov<-ggplot(dat1,aes(reorder(time,order),est,group=1)) +
  geom_hline(yintercept = 0,linetype="dashed") +
  #  geom_line() +
  geom_point() +
  geom_errorbar(aes(ymin=conf.lo,
                    ymax=conf.hi),width=0) +
  theme_r() +
  scale_y_continuous(limits=c(-1.3,1)) +
  labs(y="Average treatment effect",
       x=NULL,
       title="Government exit")  
### Matched sets
art_match<-list(dat_art$matched.sets)
art_match<-data.frame(id_year = objects(art_match[[1]]))
art_match<-art_match %>% mutate(id=gsub("\\..*","",id_year),
                                year=gsub("^.*\\.","",id_year))
art_match <- merge(art_match,cnames,by="id",all.x=T)
art_match <- art_match %>%
  mutate(cy = paste0(country," (",year,")")) %>%
  dplyr::select(cy)
paste(art_match$cy,collapse = " ")
### Combine 
f_gov<-ggarrange(f1_gov,f2_gov,widths = c(1,0.7))

# PARLIAMENT 
library(PanelMatch)
cl$id <- as.integer(as.numeric(as.factor(cl$country)))
cl$time <- as.integer(cl$year)
method <- "ps.weight"
## Parliament start
cl2 <- PanelData(panel.data = cl,
                 unit.id = "id",
                 time.id = "time",
                 treatment = "PRR_parl_dummy",
                 outcome = "Satis")
p_parl_att <- PanelMatch(cl2,
                         lag = 5,
                         refinement.method = method, 
                         match.missing = TRUE,
                         covs.formula = ~ vdem_corr + wdi_gdpgr + unemployment + wpp_netmig + I(lag(Satis, 1:4)),
                         size.match = 5, qoi = "att", 
                         lead = 0:10, forbid.treatment.reversal = FALSE,
                         use.diagonal.variance.matrix = TRUE)
dat_att<-PanelEstimate(cl2, sets = p_parl_att,
                       number.iterations = 5000,
                       confidence.level = .95)
dat1<-data.frame(est=dat_att$estimate,
                 se=dat_att$standard.error)
dat1$time <- rownames(dat1)
dat1$conf.hi <- dat1$est+dat1$se*1.96
dat1$conf.lo <- dat1$est-dat1$se*1.96
dat1$order <- seq(0,10,1)
f1_parl<-ggplot(dat1,aes(reorder(time,order),est,group=1)) +
  geom_hline(yintercept = 0,linetype="dashed") +
  geom_point() +
  geom_errorbar(aes(ymin=conf.lo,
                    ymax=conf.hi),width=0) +
  theme_r() +
  scale_y_continuous(limits=c(-0.5,1)) +  
  scale_x_discrete(breaks=c("t+0","t+2","t+4","t+6","t+8","t+10")) +  
  labs(y=NULL,
       x=NULL,
       title="Parliament entry")
### Tables for matched sets
cnames <- cl %>% 
  group_by(country) %>%
  dplyr::select(id,country) %>%
  slice(1)
### Entry
att_match<-list(dat_att$matched.sets)
att_match<-data.frame(id_year = objects(att_match[[1]]))
att_match<-att_match %>% mutate(id=gsub("\\..*","",id_year),
                                year=gsub("^.*\\.","",id_year))
att_match <- merge(att_match,cnames,by="id",all.x=T)
att_match <- att_match %>%
  mutate(cy = paste0(country," (",year,")")) %>%
  dplyr::select(cy)
paste(att_match$cy,collapse = " ")

## Parliament end
cl2 <- PanelData(panel.data = cl,
                 unit.id = "id",
                 time.id = "time",
                 treatment = "PRR_parl_dummy",
                 outcome = "Satis")
p_parl_art <- PanelMatch(cl2,
                         lag = 5,
                         refinement.method = method, 
                         match.missing = TRUE,
                         covs.formula = ~ vdem_corr + wdi_gdpgr + unemployment + wpp_netmig + I(lag(Satis, 1:2)),
                         size.match = 5, qoi = "art", 
                         lead = 0:5, forbid.treatment.reversal = FALSE,
                         use.diagonal.variance.matrix = TRUE)
dat_art<-PanelEstimate(cl2, sets = p_parl_art,
                       number.iterations = 5000,
                       confidence.level = .95)
dat1<-data.frame(est=dat_art$estimate,
                 se=dat_art$standard.error)
dat1$time <- rownames(dat1)
dat1$conf.hi <- dat1$est+dat1$se*1.96
dat1$conf.lo <- dat1$est-dat1$se*1.96
dat1$order <- seq(0,5,1)
f2_parl<-ggplot(dat1,aes(reorder(time,order),est,group=1)) +
  geom_hline(yintercept = 0,linetype="dashed") +
  #  geom_line() +
  geom_point() +
  geom_errorbar(aes(ymin=conf.lo,
                    ymax=conf.hi),width=0) +
  theme_r() +
  scale_y_continuous(limits=c(-1.3,1)) +
  labs(y=NULL,
       x=NULL,
       title="Parliament exit")  
### Matched sets
art_match<-list(dat_art$matched.sets)
art_match<-data.frame(id_year = objects(art_match[[1]]))
art_match<-art_match %>% mutate(id=gsub("\\..*","",id_year),
                                year=gsub("^.*\\.","",id_year))
art_match <- merge(art_match,cnames,by="id",all.x=T)
art_match <- art_match %>%
  mutate(cy = paste0(country," (",year,")")) %>%
  dplyr::select(cy)
paste(art_match$cy,collapse = " ")
# FIGURES
## Entry
fig_entry<-ggarrange(f1_gov,f1_parl,widths = c(1,0.95))
ggsave(fig_entry,
       file="PRR ENTRY, panel matching, claasen data.jpeg",
       width=5,
       height=2.8,
       dpi=1000)
## Exit
fig_exit<-ggarrange(f2_gov,f2_parl,widths = c(1,0.95))
ggsave(fig_exit,
       file="PRR EXIT, panel matching, claasen data.jpeg",
       width=5,
       height=2.8,
       dpi=1000)



########## ESS analysis ########## 

# Data
d <- read_dta("ESS 1-11 partyfam winlose final.dta")
d_old <- d 
#d <- d_old
d$code<-countrycode(d$cntry,"iso2c","iso3c")
d <- d %>% 
  mutate(year = coalesce(inwyye,
                         inwyr,
                         as.numeric(substr(d$inwde,1,4))))
d$satisfy <- d$stfdem
d$fam <- as.character(sjlabelled::as_label(d$partyfam))
d$fam_prr <- ifelse(d$fam=="Nationalist","PRR","Others")
d$fam[d$fam=="Nationalist"] <- "PRR"
## Demographic 
d$gender <- sjlabelled::as_label(d$gndr)
d$edu <- d$eisced
d$edu[d$edu==55] <- NA
d$age<-d$agea

## Satisfy standardized
d$satisfy_sd <- (d$satisfy-mean(d$satisfy,na.rm=T))/sd(d$satisfy,na.rm=T) 

## Merge with country-level data
d1 <- merge(d,prr,by=c("code","year"),all.x = T)
## Country_round
d1$country_round <- paste0(d1$country,"_", d1$essround)

# Filter only treated cases with data on SWD
treated_gov <- d1 %>% 
  group_by(country,PRR_government_support_lag1) %>%
  dplyr::summarise(mean=mean(satisfy_sd,na.rm=T)) %>%
  na.omit() %>%
  spread(PRR_government_support_lag1,mean) %>%
  na.omit() %>%
  dplyr::select(c=country)
treated_parl <- d1 %>% 
  group_by(country,PRR_parl_dummy_lag1) %>%
  dplyr::summarise(mean=mean(satisfy_sd,na.rm=T)) %>%
  na.omit() %>%
  spread(PRR_parl_dummy_lag1,mean) %>%
  na.omit() %>%
  dplyr::select(c=country)
d1$winner <- d1$winlose
d1_treated <- d1 %>% 
  filter(country %in% treated_gov$c)
d1_treated2 <- d1 %>% 
  filter(country %in% treated_parl$c) 


# ESS Figure (Controlling for winner loser)
d1_treated_nna <- d1_treated %>% 
  dplyr::select(satisfy_sd,
                PRR_government_support_lag1,
                fam,
                year,
                gender,
                age,
                edu,
                code,
                country_round,
                pspwght,
                winner) %>%
  mutate(fam=as.character(fam)) %>%
  na.omit()

## BEFORE AND AFTER CONTROLS FOR WINNER LOSER
### Before 
m1<-lm_robust(satisfy_sd~PRR_government_support_lag1+code+year+gender+age+edu,
              data=d1_treated_nna,
              weights = d1_treated_nna$pspwght,
              clusters = d1_treated_nna$country_round, 
              se_type = "stata")
m2<-lm_robust(satisfy_sd~PRR_government_support_lag1*fam+code+year+gender+age+edu,
              data=d1_treated_nna,
              weights = d1_treated_nna$pspwght,
              clusters = d1_treated_nna$country_round, 
              se_type = "stata")

res<-cond_effect(m2,x="PRR_government_support_lag1",w="fam")
res$coef <- res$`x's Effect`
res$se<-res$`Std. Error`
res_all <- data.frame(fam="Overall effect",
                      coef=summary(m1)$coef[2,1],
                      se=summary(m1)$coef[2,2])
res <- rbind.fill(res_all,res)
res$conf.hi <- res$coef+res$se*1.96
res$conf.lo <- res$coef-res$se*1.96
res$group <- " Baseline controls"
res_before <- res
### After 
m1<-lm_robust(satisfy_sd~PRR_government_support_lag1+code+year+gender+age+edu+winner,
              data=d1_treated_nna,
              weights = d1_treated_nna$pspwght,
              clusters = d1_treated_nna$country_round, 
              se_type = "stata")
m2<-lm_robust(satisfy_sd~PRR_government_support_lag1*fam+code+year+gender+age+edu+winner,
              data=d1_treated_nna,
              weights = d1_treated_nna$pspwght,
              clusters = d1_treated_nna$country_round, 
              se_type = "stata")
res<-cond_effect(m2,x="PRR_government_support_lag1",w="fam")
res$coef <- res$`x's Effect`
res$se<-res$`Std. Error`
res_all <- data.frame(fam="Overall effect",
                      coef=summary(m1)$coef[2,1],
                      se=summary(m1)$coef[2,2])
res <- rbind.fill(res_all,res)
res$conf.hi <- res$coef+res$se*1.96
res$conf.lo <- res$coef-res$se*1.96
res_after <- res
res_after$group <- "+ winner/loser status"
## Combine and Change order
res <- rbind(res_before,res_after)
res$order <- rank(res$coef)
res$order[res$fam=="Overall effect"] <- 99
res$fam<-reorder(res$fam,res$order)
res <- res %>% filter(!c(fam=="PRR" & group=="+ winner/loser status"))
## Plot
f1<-ggplot(res,aes(fam,coef,color=fam,shape=group)) +
  geom_hline(yintercept = 0) +
  geom_errorbar(aes(ymin=conf.lo,
                    ymax=conf.hi),
                width=0.2,
                position = position_dodge(width=-0.4)) +
  geom_point(position = position_dodge(width=-0.4),fill="white") +  
  theme_r() +  
  theme(legend.position = c(0.75,0.15)) +
  scale_color_manual(values=c(pal,"black")) +
  scale_shape_manual(values=c(16,13)) +
  labs(x=NULL,
       y="Effect of PRR government inclusion\non satisfaction with democracy") +
  coord_flip() +
  labs(shape=NULL) +
  guides(color = "none",
         shape=guide_legend(nrow=2))
## Export
ggsave(f1,
       file="PRR government inclusion, ESS, effect by party (control for winner loser).jpeg",
       width=5,
       height=4,
       dpi=1000)








############################################################
######################### APPENDIX #########################
############################################################


# APX Table desc stats 
cl_sub <- cl %>% dplyr::select(Satis,
                               PRR_government_support,
                               PRR_parl_dummy,
                               vdem_corr,
                               wdi_gdpgr,
                               unemployment,
                               wpp_netmig)
datasummary_skim(cl_sub,
                 output = "APX desc table.docx")




###################### Synthetic Control ######################
library(gsynth)
library(tidyr)
library(plyr)
system.time(out <- gsynth(Satis ~ 
                            PRR_government_support +
                            vdem_corr + 
                            unemployment +
                            wdi_gdpgr + 
                            wpp_netmig,
                          data = cl,
                          index = c("country","year"), se = TRUE,force="two-way",
                          min.T0 = 10,
                          CV = F,
                          na.rm=T,
                          inference = "parametric", nboots = 1000, 
                          parallel = TRUE))
p<-plot(out)
data<-p$data
stats<-data.frame(round(t(data.frame(out$est.avg)),3))
f1<-ggarrange(ggplot(data, aes(x=time)) +
                geom_vline(xintercept=0,alpha=0.3,linetype="dashed") +
                geom_hline(yintercept=0,) +
                geom_line(aes(y=ATT),color="black",linetype="solid",size=0.7) +
                geom_ribbon(aes(ymin=CI.lower,ymax=CI.upper),alpha=0.2) +
                scale_x_continuous(limits=c(-20,20)) +
                scale_y_continuous(limits=c(-0.75,2.2)) +
                theme_classic() +
                labs(x="Years relative to treatment (Government Inclusion)",
                     y="Effect on satisfaction with democracy") +
                theme(panel.grid = element_blank(),
                      strip.background = element_blank(),
                      text = element_text(color="black"),
                      axis.text = element_text(color="black"))) +
  geom_text(aes(x=0.15,y=0.9, hjust=0,
                label="Average Treatment Effect"),
            fontface="bold",
            size=4) +
  geom_text(aes(x=0.15,y=0.8, hjust=0,
                label=paste0(
                  "\nEstimate: ", stats[1,],
                  "\nStd. Err: ", stats[2,],
                  "\np-value: <0.001")),
            size=4)
ggsave(f1,
       file="APX PRR, synthetic control, claasen data.jpeg",
       width=4,
       height=3.2,
       dpi=1000)

##### Synthetic: Parliament
system.time(out <- gsynth(Satis ~ 
                            PRR_parl_dummy +
                            vdem_corr + 
                            unemployment +
                            wdi_gdpgr + 
                            wpp_netmig, 
                          data = cl,
                          index = c("country","year"), se = TRUE,force="two-way",
                          min.T0 = 7,
                          CV = F,
                          na.rm=T,
                          inference = "parametric", nboots = 1000, 
                          parallel = TRUE))
p<-plot(out)
data<-p$data
stats<-data.frame(round(t(data.frame(out$est.avg)),3))

f1<-ggarrange(ggplot(data, aes(x=time)) +
                geom_vline(xintercept=0,alpha=0.3,linetype="dashed") +
                geom_hline(yintercept=0,) +
                geom_line(aes(y=ATT),color="black",linetype="solid",size=0.7) +
                geom_ribbon(aes(ymin=CI.lower,ymax=CI.upper),alpha=0.2) +
                scale_x_continuous(limits=c(-20,20)) +
                #  scale_y_continuous(limits=c(-2,1.5)) +
                theme_classic() +
                labs(x="Years relative to treatment (Parliament Inclusion)",
                     y="Effect on satisfaction with democracy") +
                theme(panel.grid = element_blank(),
                      strip.background = element_blank(),
                      text = element_text(color="black"),
                      axis.text = element_text(color="black"))) +
  geom_text(aes(x=0.15,y=0.9, hjust=0,
                label="Average Treatment Effect"),
            fontface="bold",
            size=4) +
  geom_text(aes(x=0.15,y=0.8, hjust=0,
                label=paste0(
                  "\nEstimate: ", stats[1,],
                  "\nStd. Err: ", stats[2,],
                  "\np-value: ", stats[5,])),
            size=4)
ggsave(f1,
       file="APX PRR, synthetic control PARLIAMENT, claasen data.jpeg",
       width=4,
       height=3.2,
       dpi=1000)


# APX FIGURE: Covariate balance 
covs <- c("vdem_corr","wdi_gdpgr","unemployment","wpp_netmig","Satis")
## Government entry
bal<-get_covariate_balance(p_govsup_att,
                           panel.data = cl1,
                           covariates = covs)
n<-(as.numeric(nrow(data.frame(bal)))-1)*-1
dat <- rbind(data.frame(get_unrefined_balance(bal)[[1]],
                        group=" Before refinement",
                        t=seq(n,0,1)),
             data.frame(bal[[1]],
                        group="After refinement",
                        t=seq(n,0,1)))

dat1 <- dat %>%
  gather(var,val,att.vdem_corr:att.Satis) %>%
  mutate(var = gsub("^.*\\.","",var),
         treatment="Government entry")

## Government exit
bal<-get_covariate_balance(p_govsup_art,
                           panel.data = cl1,
                           covariates = covs)
n<-(as.numeric(nrow(data.frame(bal)))-1)*-1
dat <- rbind(data.frame(get_unrefined_balance(bal)[[1]],
                        group=" Before refinement",
                        t=seq(n,0,1)),
             data.frame(bal[[1]],
                        group="After refinement",
                        t=seq(n,0,1)))
dat2 <- dat %>%
  gather(var,val,art.vdem_corr:art.Satis) %>%
  mutate(var = gsub("^.*\\.","",var),
         treatment="Government exit")

## Parliament entry
bal<-get_covariate_balance(p_parl_att,
                           panel.data = cl2,
                           covariates = covs)
n<-(as.numeric(nrow(data.frame(bal)))-1)*-1
dat <- rbind(data.frame(get_unrefined_balance(bal)[[1]],
                        group=" Before refinement",
                        t=seq(n,0,1)),
             data.frame(bal[[1]],
                        group="After refinement",
                        t=seq(n,0,1)))
dat3 <- dat %>%
  gather(var,val,att.vdem_corr:att.Satis) %>%
  mutate(var = gsub("^.*\\.","",var),
         treatment="Parliament entry")

## Parliament exit
bal<-get_covariate_balance(p_parl_art,
                           panel.data = cl2,
                           covariates = covs)
n<-(as.numeric(nrow(data.frame(bal)))-1)*-1
dat <- rbind(data.frame(get_unrefined_balance(bal)[[1]],
                        group=" Before refinement",
                        t=seq(n,0,1)),
             data.frame(bal[[1]],
                        group="After refinement",
                        t=seq(n,0,1)))
dat4 <- dat %>%
  gather(var,val,art.vdem_corr:art.Satis) %>%
  mutate(var = gsub("^.*\\.","",var),
         treatment="Parliament exit")

# Combine and plot 
dat <- rbind.fill(dat1,dat2,dat3,dat4)

f1<-ggplot(dat,aes(t,val,color=var)) +
  geom_line() +
  facet_grid(treatment~group) +
  scale_y_continuous(limits = c(-2,2)) +
  theme_r() +
  labs(y="Average difference between treated and control\n(in standard deviations)") +
  scale_color_npg()

ggsave(f1,file="APX prr, covariate balance.jpeg",
       width=6,height=6,
       dpi=1000)



# APX FIGURE: Treatment history
f1<-ggplot(cl,aes(year,reorder(country,PRR_government_support),fill=factor(PRR_government_support))) +
  geom_tile(color="white") +
  theme_r() +
  theme(legend.position = "none") +
  scale_fill_manual(values=c("grey","darkblue")) +
  labs(x=NULL,y=NULL)
ggsave(f1,
       file="APX PRR, government treatment history.jpeg",
       width=6,height=5,
       dpi=1000)
f1<-ggplot(cl,aes(year,reorder(country,PRR_parl_dummy),fill=factor(PRR_parl_dummy))) +
  geom_tile(color="white") +
  theme_r() +
  theme(legend.position = "none") +
  scale_fill_manual(values=c("grey","darkblue")) +
  labs(x=NULL,y=NULL)
ggsave(f1,
       file="APX PRR, parliament treatment history.jpeg",
       width=6,height=5,
       dpi=1000)



########## Regression analysis ###########

# Split sample 
cl_subset <- cl %>% filter(year>=2000)
m1 <- lm(Satis~PRR_government_support_lag1 +
           vdem_corr + 
           unemployment +
           wdi_gdpgr + 
           wpp_netmig+
           country,
         data=cl_subset)
cl_subset <- cl %>% filter(year<2000)
m2 <- lm(Satis~PRR_government_support_lag1 +
           vdem_corr + 
           unemployment +
           wdi_gdpgr + 
           wpp_netmig+
           country,
         data=cl_subset)
modelsummary(list("2000-2020"=m1,"1973-1999"=m2),
             output = "gt",
             stars = T, 
             coef_map = c("PRR_government_support_lag1"="PRR government inclusion",
                          "PRR_parl_dummy_lag1" = "PRR parliament inclusion",
                          "vdem_corr" = "Corruption (0-1)",
                          "unemployment" = "Unemployment (%)",
                          "wdi_gdpgr" = "GDP growth (%)",
                          "wpp_netmig" = "Net immigration (per 1000)",
                          "(Intercept)" = "Intercept"),
             gof_map = c("nobs","r.squared"),
) %>%
  tab_spanner(label="Dependent variable: Satisfaction with democracy",
              columns=c(2:3)) %>%  
  gtsave(filename="APX prr reg table 2000.docx")




########## ESS analysis ###########

# Effect of PRR parliament inclusion by party
d1_treated_nna <- d1_treated2 %>% 
  dplyr::select(satisfy_sd,
                PRR_parl_dummy_lag1,
                fam,
                year,
                gender,
                age,
                edu,
                code,
                country_round,
                pspwght) %>%
  na.omit()
m1<-lm_robust(satisfy_sd~PRR_parl_dummy_lag1+code+year+gender+age+edu,
              data=d1_treated_nna,
              weights = d1_treated_nna$pspwght,
              clusters = d1_treated_nna$country_round, 
              se_type = "stata")
m2<-lm_robust(satisfy_sd~PRR_parl_dummy_lag1*fam+code+year+gender+age+edu,
              data=d1_treated_nna,
              weights = d1_treated_nna$pspwght,
              clusters = d1_treated_nna$country_round, 
              se_type = "stata")
# Effect by party family
res<-cond_effect(m2,x="PRR_parl_dummy_lag1",w="fam")
res$coef <- res$`x's Effect`
res$se<-res$`Std. Error`
# Add overall effect 
res_all <- data.frame(fam="Overall effect",
                      coef=summary(m1)$coef[2,1],
                      se=summary(m1)$coef[2,2])
res <- rbind.fill(res_all,res)
# Calculate CIs 
res$conf.hi <- res$coef+res$se*1.96
res$conf.lo <- res$coef-res$se*1.96
# Change order
res$order <- rank(res$coef)
res$order[res$fam=="Overall effect"] <- 12
res$fam<-reorder(res$fam,res$order)
# Plot
f1<-ggplot(res,aes(fam,coef,fill=fam)) +
  geom_bar(stat="identity",color="black") +
  geom_errorbar(aes(ymin=conf.lo,
                    ymax=conf.hi),
                width=0.2) +
  theme_r() +  
  theme(legend.position = "none") +
  scale_fill_manual(values=c(pal,"white")) +
  labs(x=NULL,
       y="Effect of PRR parliament inclusion\non satisfaction with democracy") +
  coord_flip() 
ggsave(f1,
       file="APX PRR parliament inclusion, ESS, effect by party.jpeg",
       width=4,
       height=3.2,
       dpi=1000)


# Effect of PRR government inclusion by country
d1_treated_nna <- d1_treated %>% 
  dplyr::select(satisfy_sd,
                PRR_government_support_lag1,
                fam,
                year,
                gender,
                age,
                edu,
                country,
                country_round,
                pspwght) %>%
  na.omit()


m1<-lm_robust(satisfy_sd~PRR_government_support_lag1+country+year+gender+age+edu,
              data=d1_treated_nna,
              weights = d1_treated_nna$pspwght,
              clusters = d1_treated_nna$country_round, 
              se_type = "stata")
m2<-lm_robust(satisfy_sd~PRR_government_support_lag1*country+country+year+gender+age+edu,
              data=d1_treated_nna,
              weights = d1_treated_nna$pspwght,
              clusters = d1_treated_nna$country_round, 
              se_type = "stata")
# Effect by party family
res<-cond_effect(m2,x="PRR_government_support_lag1",w="country")
res$coef <- res$`x's Effect`
res$se<-res$`Std. Error`
# Add overall effect 
res_all <- data.frame(country="Overall effect",
                      coef=summary(m1)$coef[2,1],
                      se=summary(m1)$coef[2,2])
res <- rbind.fill(res_all,res)
# Calculate CIs 
res$conf.hi <- res$coef+res$se*1.96
res$conf.lo <- res$coef-res$se*1.96
# Change order
res$order <- rank(res$coef)
res$order[res$country=="Overall effect"] <- 99
res$country<-reorder(res$country,res$order)
# Plot
f1<-ggplot(res,aes(country,coef,fill=country)) +
  geom_bar(stat="identity",color="black") +
  geom_errorbar(aes(ymin=conf.lo,
                    ymax=conf.hi),
                width=0.2) +
  theme_r() +  
  theme(legend.position = "none") +
  scale_fill_manual(values=c(pal,pal[c(1,2,3)],"white")) +
  labs(x=NULL,
       y="Effect of PRR government inclusion\non satisfaction with democracy") +
  coord_flip() 

ggsave(f1,
       file="APX PRR inclusion government, ESS, effect by country.jpeg",
       width=4,
       height=3.2,
       dpi=1000)


# Effect of PRR parliament inclusion by country
d1_treated_nna <- d1_treated2 %>% 
  dplyr::select(satisfy_sd,
                PRR_parl_dummy_lag1,
                fam,
                year,
                gender,
                age,
                edu,
                code,
                country,
                country_round,
                pspwght) %>%
  na.omit()

m1<-lm_robust(satisfy_sd~PRR_parl_dummy_lag1+country+year+gender+age+edu,
              data=d1_treated_nna,
              weights = d1_treated_nna$pspwght,
              clusters = d1_treated_nna$country_round, 
              se_type = "stata")
m2<-lm_robust(satisfy_sd~PRR_parl_dummy_lag1*country+year+gender+age+edu,
              data=d1_treated_nna,
              weights = d1_treated_nna$pspwght,
              clusters = d1_treated_nna$country_round, 
              se_type = "stata")
# Effect by party family
res<-cond_effect(m2,x="PRR_parl_dummy_lag1",w="country")
res$coef <- res$`x's Effect`
res$se<-res$`Std. Error`
# Add overall effect 
res_all <- data.frame(country="Overall effect",
                      coef=summary(m1)$coef[2,1],
                      se=summary(m1)$coef[2,2])
res <- rbind.fill(res_all,res)
# Calculate CIs 
res$conf.hi <- res$coef+res$se*1.96
res$conf.lo <- res$coef-res$se*1.96
# Change order
res$order <- rank(res$coef)
res$order[res$country=="Overall effect"] <- 99
res$country<-reorder(res$country,res$order)
# Plot
f1<-ggplot(res,aes(country,coef,fill=country)) +
  geom_bar(stat="identity",color="black") +
  geom_errorbar(aes(ymin=conf.lo,
                    ymax=conf.hi),
                width=0.2) +
  theme_r() +  
  theme(legend.position = "none") +
  scale_fill_manual(values=c(pal,pal[c(1,2)],"white")) +
  labs(x=NULL,
       y="Effect of PRR parliament inclusion\non satisfaction with democracy") +
  coord_flip() 


ggsave(f1,
       file="APX PRR inclusion parliament, ESS, effect by country.jpeg",
       width=4,
       height=3.2,
       dpi=1000)



# Regression results for ESS Figure in manuscript
# Government
d1_treated_nna <- d1_treated %>% 
  dplyr::select(satisfy_sd,
                PRR_government_support_lag1,
                fam,
                year,
                gender,
                age,
                edu,
                code,
                country_round,
                pspwght,
                winner) %>%
  na.omit()

m1_gov<-lm(satisfy_sd~PRR_government_support_lag1+code+year+gender+age+edu+winner,
           data=d1_treated_nna,weights=pspwght)
m2_gov<-lm(satisfy_sd~PRR_government_support_lag1*fam+code+year+gender+age+edu+winner,
           data=d1_treated_nna,weights=pspwght)

library(kableExtra)
library(gt)
# Table
msummary(list(m1_gov,m2_gov),
         vcov=~country_round,
         stars=T,
         output = "table with winner.docx",
         gof_map = c("nobs","r.squared"),
         coef_rename = c("PRR_government_support_lag1"="Inclusion (gov/parl)",
                         "PRR_parl_dummy_lag1"="Inclusion (gov/parl)",
                         "year" = "Year",
                         "age" = "Age",
                         "genderFemale" = "Female",
                         "edu" = "Education",
                         "winner" = "Electoral winner (0,1)"),
         coef_omit = c("code"))

tidy_with_ref <- function(model, ...) {
  tidy(model, ...) %>%
    tidy_add_reference_rows()
}

modelsummary(m1_gov,tidy=tidy_with_ref)



########## Jackknifing summary ##########

### CL REGRESSION 
# Initialize an empty dataframe to store results
jackknife_results <- data.frame(
  country_removed = character(),
  coef = numeric(),
  stringsAsFactors = FALSE)
countries <- c(unique(cl$country),"Overall")
for (ctry in countries) {
  # Subset data, excluding current country
  cl_subset <- subset(cl, country != ctry)
  # Run analysis
  m <- lm(Satis~PRR_government_support_lag1 +
            vdem_corr + 
            unemployment +
            wdi_gdpgr + 
            wpp_netmig+
            country,
          data=cl_subset)
  # Store coef 
  coef=summary(m)$coef[2,1]
  # Append result to dataframe
  jackknife_results <- rbind(jackknife_results, 
                             data.frame(country_removed = ctry,
                                        coef = coef))
}
jack_reg <- jackknife_results


#### PANEL MATCHING
method <- "ps.weight"
# Initialize an empty dataframe to store results
jackknife_results <- data.frame(
  country_removed = character(),
  mean_estimate = numeric(),
  stringsAsFactors = FALSE
)
countries <- c(unique(cl$country),"Overall")
# Loop over unique countries
for (ctry in countries) {
  
  # Subset data, excluding current country
  cl_subset <- subset(cl, country != ctry)
  
  # Create PanelData object
  cl1 <- PanelData(panel.data = cl_subset,
                   unit.id = "id",
                   time.id = "time",
                   treatment = "PRR_government_support",
                   outcome = "Satis")
  
  # Run PanelMatch
  p_govsup_att <- PanelMatch(cl1,
                             lag = 5,
                             refinement.method = method, 
                             match.missing = TRUE,
                             covs.formula = ~ vdem_corr + wdi_gdpgr + unemployment + wpp_netmig + I(lag(Satis, 1:4)),
                             size.match = 5, qoi = "att", 
                             lead = 0:10, forbid.treatment.reversal = FALSE,
                             use.diagonal.variance.matrix = TRUE)
  
  # Estimate ATT
  dat_att <- PanelEstimate(cl1, sets = p_govsup_att,
                           number.iterations = 10,
                           confidence.level = .95)
  
  # Compute mean estimate (excluding the first element)
  mean_est <- mean(dat_att$estimate[-c(1:3)])
  
  # Append result to dataframe
  jackknife_results <- rbind(jackknife_results, 
                             data.frame(country_removed = ctry,
                                        mean_estimate = mean_est))
}
jack_panel<-jackknife_results

### ESS ANALYSIS 
# Initialize an empty dataframe to store results
jackknife_results <- data.frame(
  country_removed = character(),
  coef = numeric(),
  stringsAsFactors = FALSE)

d1_treated_nna <- d1_treated %>% 
  dplyr::select(satisfy_sd,
                PRR_government_support_lag1,
                fam,
                year,
                gender,
                age,
                edu,
                code,
                country_round,
                pspwght,
                winner) %>%
  mutate(fam=as.character(fam)) %>%
  na.omit()
countries <- c(unique(d1_treated_nna$code),"Overall")

for (ctry in countries) {
  # Subset data, excluding current country
  d_subset <- subset(d1_treated_nna, code != ctry)
  # Run analysis
  m<-lm(satisfy_sd~PRR_government_support_lag1+code+year+gender+age+edu+winner,
        data=d_subset,weights = d_subset$pspwght)
  # Store coef 
  coef=summary(m)$coef[2,1]
  # Append result to dataframe
  jackknife_results <- rbind(jackknife_results, 
                             data.frame(country_removed = ctry,
                                        coef = coef))
}
jack_ess <- jackknife_results

dat<-data.frame(analysis = c("Regression with aggregate data (coefficient)",
                             "Panel matching (average treatment effect, t3-t10)",
                             "Regression with individual data (coefficient)"),
                estimate=c(jack_reg$coef[jack_reg$country_removed=="Overall"],
                           jack_panel$mean_estimate[jack_panel$country_removed=="Overall"],
                           jack_ess$coef[jack_ess$country_removed=="Overall"]),
                min=c(min(jack_reg$coef),
                      min(jack_panel$mean_estimate),
                      min(jack_ess$coef)),
                max=c(max(jack_reg$coef),
                      max(jack_panel$mean_estimate),
                      max(jack_ess$coef)))
dat[,2:4] <- round(dat[,2:4],3)

datasummary_df(dat, output = "jackknife.docx")




########## Jackknifing (distribution of results) ##########

### --- CL REGRESSION --- ###
jackknife_results <- data.frame(
  country_removed = character(),
  coef = numeric(),
  stringsAsFactors = FALSE
)

countries <- c(unique(cl$country), "Overall")

for (ctry in countries) {
  cl_subset <- subset(cl, country != ctry)
  m <- lm(Satis ~ PRR_government_support_lag1 +
            vdem_corr + unemployment + wdi_gdpgr +
            wpp_netmig + country,
          data = cl_subset)
  
  coef <- summary(m)$coef[2,1]
  
  jackknife_results <- rbind(jackknife_results,
                             data.frame(country_removed = ctry,
                                        coef = coef))
}
jack_reg <- jackknife_results


### --- PANEL MATCHING --- ###
method <- "ps.weight"

jackknife_results <- data.frame(
  country_removed = character(),
  mean_estimate = numeric(),
  stringsAsFactors = FALSE
)

countries <- c(unique(cl$country), "Overall")

for (ctry in countries) {
  
  cl_subset <- subset(cl, country != ctry)
  
  cl1 <- PanelData(
    panel.data = cl_subset,
    unit.id = "id",
    time.id = "time",
    treatment = "PRR_government_support",
    outcome = "Satis"
  )
  
  p_govsup_att <- PanelMatch(
    cl1,
    lag = 5,
    refinement.method = method,
    match.missing = TRUE,
    covs.formula = ~ vdem_corr + wdi_gdpgr + unemployment +
      wpp_netmig + I(lag(Satis, 1:4)),
    size.match = 5,
    qoi = "att",
    lead = 0:10,
    forbid.treatment.reversal = FALSE,
    use.diagonal.variance.matrix = TRUE
  )
  
  dat_att <- PanelEstimate(
    cl1,
    sets = p_govsup_att,
    number.iterations = 10,
    confidence.level = .95
  )
  
  mean_est <- mean(dat_att$estimate[-c(1:3)])
  
  jackknife_results <- rbind(
    jackknife_results,
    data.frame(country_removed = ctry,
               mean_estimate = mean_est)
  )
}

jack_panel <- jackknife_results


### --- ESS ANALYSIS --- ###
jackknife_results <- data.frame(
  country_removed = character(),
  coef = numeric(),
  stringsAsFactors = FALSE
)

d1_treated_nna <- d1_treated %>%
  dplyr::select(satisfy_sd, PRR_government_support_lag1, fam, year,
                gender, age, edu, code, country_round,
                pspwght, winner) %>%
  mutate(fam = as.character(fam)) %>%
  na.omit()

countries <- c(unique(d1_treated_nna$code), "Overall")

for (ctry in countries) {
  
  d_subset <- subset(d1_treated_nna, code != ctry)
  
  m <- lm(satisfy_sd ~ PRR_government_support_lag1 + code + year +
            gender + age + edu + winner,
          data = d_subset,
          weights = d_subset$pspwght)
  
  coef <- summary(m)$coef[2,1]
  
  jackknife_results <- rbind(jackknife_results,
                             data.frame(country_removed = ctry,
                                        coef = coef))
}

jack_ess <- jackknife_results


### --- COMBINED OUTPUT: ALL VALUES --- ###

# Long-format dataset containing all jackknife results
all_jackknife <- rbind(
  data.frame(
    analysis = "Regression (Aggregate)",
    country_removed = jack_reg$country_removed,
    value = jack_reg$coef
  ),
  data.frame(
    analysis = "Panel Matching",
    country_removed = jack_panel$country_removed,
    value = jack_panel$mean_estimate
  ),
  data.frame(
    analysis = "Regression (Individual)",
    country_removed = jack_ess$country_removed,
    value = jack_ess$coef
  )
)

dat <- all_jackknife %>%
  dplyr::group_by(analysis) %>%
  dplyr::summarise(
    overall = value[country_removed == "Overall"],
    n_runs = n(),
    .groups = "drop"
  )

dat[,2:ncol(dat)] <- round(dat[,2:ncol(dat)], 3)

datasummary_df(dat, output = "jackknife.docx")

f1<-ggplot(all_jackknife,aes(value)) +
  geom_vline(xintercept = 0,color="white") +
  geom_histogram() +
  facet_wrap(~analysis,scales="free") +
  theme_r() + 
  theme(panel.spacing = unit(2, "lines")) +
  coord_cartesian(clip="off") +
  labs(x="Estimated effect of PRR government inclusion on satisfaction with democracy (standardized)",
       y="Number of iterations")

ggsave(f1,
       file="APX jackknife distributions.jpeg",
       width=6,height=2.5,
       dpi=1000)






# CONTROLLING FOR VOTE SHARE, GROWTH, AND THEIR INTERACTION

## Country-year regression analysis
m1 <- lm_robust(Satis~PRR_government_support_lag1+
                  code,
                data=cl,
                clusters = code,
                se_type="stata")
m2 <- lm_robust(Satis~PRR_government_support_lag1 +
                  vdem_corr + 
                  unemployment +
                  wdi_gdpgr + 
                  wpp_netmig+
                  code,
                data=cl,
                clusters = code,
                se_type="stata")
m3 <- lm_robust(Satis~PRR_government_support_lag1 +
                  vdem_corr + 
                  unemployment +
                  wdi_gdpgr + 
                  wpp_netmig+
                  PRR_vote_share+
                  code,
                data=cl,
                clusters = code,
                se_type="stata")
m4 <- lm_robust(Satis~PRR_government_support_lag1 +
                  vdem_corr + 
                  unemployment +
                  wdi_gdpgr + 
                  wpp_netmig+
                  PRR_vote_share+
                  PRR_vote_share_growth+
                  code,
                data=cl,
                clusters = code,
                se_type="stata")
m5 <- lm_robust(Satis~PRR_government_support_lag1 +
                  vdem_corr + 
                  unemployment +
                  wdi_gdpgr + 
                  wpp_netmig+
                  PRR_vote_share*PRR_vote_share_growth+
                  code,
                data=cl,
                clusters = code,
                se_type="stata")

modelsummary(list(m1,m2,m3,m4,m5),
             output = "gt",
             stars = T, 
             coef_map = c("PRR_government_support_lag1"="PRR government inclusion",
                          "PRR_parl_dummy_lag1" = "PRR parliament inclusion",
                          "vdem_corr" = "Corruption (0-1)",
                          "unemployment" = "Unemployment (%)",
                          "wdi_gdpgr" = "GDP growth (%)",
                          "wpp_netmig" = "Net immigration (per 1000)",
                          "PRR_vote_share" = "PRR vote share (%)",
                          "PRR_vote_share_growth" ="PRR vote share growth (%-points)",
                          "PRR_vote_share:PRR_vote_share_growth"="Vote share x growth",
                          "(Intercept)" = "Intercept"),
             gof_map = c("nobs","r.squared")) %>%
  tab_spanner(label="Dependent variable: Satisfaction with democracy",
              columns=c(2:6)) %>%  
  gtsave(filename="APX country-year regression (controlling for size and growth).docx")

# Individual level regression (ESS)
m1<-lm(satisfy_sd~PRR_government_support_lag1+code+year+gender+age+edu+winner,
       data=d1_treated,weights = pspwght)
m2<-lm(satisfy_sd~PRR_government_support_lag1+code+year+gender+age+edu+winner+PRR_vote_share,
       data=d1_treated,weights = pspwght)
m3<-lm(satisfy_sd~PRR_government_support_lag1+code+year+gender+age+edu+winner+PRR_vote_share_growth,
       data=d1_treated,weights = pspwght)
m4<-lm(satisfy_sd~PRR_government_support_lag1+code+year+gender+age+edu+winner+PRR_vote_share+PRR_vote_share_growth,
       data=d1_treated,weights = pspwght)
m5<-lm(satisfy_sd~PRR_government_support_lag1+code+year+gender+age+edu+winner+PRR_vote_share*PRR_vote_share_growth,
       data=d1_treated,weights = pspwght)
modelsummary(list(m1,m2,m3,m4,m5),
             vcov=~country_round,
             stars=T,
             output = "gt",
             gof_map = c("r.squared","nobs"),
             coef_rename = c("PRR_government_support_lag1"="Government inclusion",
                             "year" = "Year",
                             "age" = "Age",
                             "genderFemale" = "Female",
                             "edu" = "Education",
                             "winner" = "Winner-loser (winner)",
                             "PRR_vote_share" = "PRR vote share (%)",
                             "PRR_vote_share_growth" ="PRR vote share growth (%-points)",
                             "PRR_vote_share:PRR_vote_share_growth"="Vote share x growth",
                             "(Intercept)" = "Intercept"),
             coef_omit = "code") %>%
  tab_spanner(label="Dependent variable: Satisfaction with democracy",
              columns=c(2:6)) %>%
  gtsave(filename="APX ESS regression (controlling for size and growth).docx")

# #Panel matching 
cl$id <- as.integer(as.numeric(as.factor(cl$country)))
cl$time <- as.integer(cl$year)
method <- "ps.weight"
### Government support start
cl1 <- PanelData(panel.data = cl,
                 unit.id = "id",
                 time.id = "time",
                 treatment = "PRR_government_support",
                 outcome = "Satis")
p_govsup_att <- PanelMatch(cl1,
                           lag = 5,
                           refinement.method = method, 
                           match.missing = TRUE,
                           covs.formula = ~ vdem_corr + wdi_gdpgr + unemployment + wpp_netmig + PRR_vote_share*PRR_vote_share_growth + I(lag(Satis, 1:4)),
                           size.match = 5, qoi = "att", 
                           lead = 0:10, forbid.treatment.reversal = FALSE,
                           use.diagonal.variance.matrix = TRUE)
dat_att<-PanelEstimate(cl1, sets = p_govsup_att,
                       number.iterations = 5000,
                       confidence.level = .95)
dat1<-data.frame(est=dat_att$estimate,
                 se=dat_att$standard.error)
dat1$time <- rownames(dat1)
dat1$conf.hi <- dat1$est+dat1$se*1.96
dat1$conf.lo <- dat1$est-dat1$se*1.96
dat1$order <- seq(0,10,1)
f1_gov<-ggplot(dat1,aes(reorder(time,order),est,group=1)) +
  geom_hline(yintercept = 0,linetype="dashed") +
  # geom_line() +
  geom_point() +
  geom_errorbar(aes(ymin=conf.lo,
                    ymax=conf.hi),width=0) +
  theme_r() +
  scale_y_continuous(limits=c(-0.5,1)) +  
  scale_x_discrete(breaks=c("t+0","t+2","t+4","t+6","t+8","t+10")) +
  labs(y="Average treatment effect",
       x=NULL,
       title="Government entry")

ggsave(f1_gov,
       file="APX panel matching (controlling for size and growth).jpeg",
       width=4,height=4,
       dpi=1000)



# APX Table desc ESS  
d1_treated_nna <- d1_treated %>% 
  dplyr::select(satisfy_sd,
                PRR_government_support_lag1,
                fam,
                gender,
                age,
                country,
                edu,
                pspwght,
                winner) %>%
datasummary_skim(d1_treated_nna,
                 output = "APX desc table ESS.docx")



# ESS winner-loser interaction (three-way) 
d1_treated_nna <- d1_treated %>% 
  dplyr::select(satisfy_sd,
                PRR_government_support_lag1,
                fam,
                year,
                gender,
                age,
                edu,
                code,
                country_round,
                pspwght,
                winner) %>%
  na.omit()
d1_treated_nna$fam_winner <- paste0(d1_treated_nna$fam,"_",d1_treated_nna$winner)

m1<-lm_robust(satisfy_sd~PRR_government_support_lag1+code+year+gender+age+edu,
              data=d1_treated_nna,
              weights = d1_treated_nna$pspwght,
              clusters = d1_treated_nna$country_round, 
              se_type = "stata")
m2<-lm_robust(satisfy_sd~PRR_government_support_lag1*fam_winner+code+year+gender+age+edu,
              data=d1_treated_nna,
              weights = d1_treated_nna$pspwght,
              clusters = d1_treated_nna$country_round, 
              se_type = "stata")
# Effect by party family
res<-cond_effect(m2,x="PRR_government_support_lag1",w="fam_winner")
res$coef <- res$`x's Effect`
res$se<-res$`Std. Error`
res$fam <- sub("_.*", "", res$fam_winner)
res$winner <- sub(".*_", "", res$fam_winner)

# Add overall effect 
res_all <- data.frame(fam="Overall effect",
                      winner=1,
                      coef=summary(m1)$coef[2,1],
                      se=summary(m1)$coef[2,2])
#res <- rbind.fill(res_all,res)
# Calculate CIs 
res$conf.hi <- res$coef+res$se*1.96
res$conf.lo <- res$coef-res$se*1.96
# Change order
res$order <- rank(res$coef)
#res$order[res$fam=="Overall effect"] <- 12
res$fam<-reorder(res$fam,res$order)

# Plot
f1<-ggplot(res%>%filter(fam!="PRR"),
           aes(fam,coef,color=winner,shape=winner)) +
  geom_point() +
  geom_errorbar(aes(ymin=conf.lo,
                    ymax=conf.hi),
                width=0.2) +
  theme_r() +  
  scale_fill_manual(values=c(pal,"white")) +
  labs(x=NULL,
       y="Effect of PRR government inclusion\non satisfaction with democracy") +
  coord_flip() 

ggsave(f1,
       file="APX family interacted with winner loser.jpeg",
       width=5,
       height=4,
       dpi=1000)


